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Abstract 


Decays of superheavy X particles with mass Mx ~ 10^^ —10^® GeV have been proposed 
as origin of the observed ultra high energy cosmic rays (UHECR). We describe in detail the 
physics involved in the different steps of the decay of such a particle. In particular, we give 
for the first time the complete set of splitting functions needed to model a parton shower 
in the minimal supersymmetric extension of the Standard Model (MSSM). We present our 
results in the form of fragmentation functions of any (s)particle of the MSSM into any hnal 
stable particle (proton, photon, electron, three types of neutrino, lightest superparticle 
LSP) at a virtuality Q = Mx, over a scaled energy range x = 2E/Mx G [10“^^,1]. 
Extending the coverage to such small fractional energies is necessary since the energy 
region around 10^® eV and below could be of considerable interest in testing this kind of 
model for generating UHECR. We explicitly demonstrate that our treatment conserves 
energy, and discuss the dependence of the final result on SUSY parameters. We also show 
that our results are essentially independent of the necessary extrapolation of the input 
fragmentation functions, which are known only for x > 0.1, towards small x. Finally, we 
added a new treatment of the color coherence effects at very small x, using the analytic 
“MLLA” solution. Our computer code will soon be made available. 
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1 Introduction 


In the last few decades ultrahigh energy cosmic rays (UHECR), with energy above the so- 
called GZK cnt-off P], have been observed [2]. The classical “bottom-np” explanation for the 
acceleration of CR particles exploits the electromagnetic helds that are likely to be present 
in objects like gamma ray bursters [Sj, “hot spots” of radio-galaxies |1] or near super-massive 
black holes in dormant quasars [H]. However, it is difficult to hnd objects capable of accelerating 
protons to energies above 10^° eV, partly because the product of held strength and spatial 
extension of the held does not seem to be sufficiently large, and partly because the accelerated 
particles can loose a fair fraction of their energy in synchrotron radiation. As an attractive 
alternative, “top-down” models have been proposed. Here one postulates the existence of 
“new physics” at a very high energy scale, i.e. the existence of super-heavy particles of masses 
greater than 10^^ GeV, which could decay and hence produce the observed UHEGR; this idea 
has been reviewed e.g. in The existence of such a very high energy scale strongly suggests 

the existence of superparticles with masses not much above 1 TeV in order to guarantee the 
perturbative stability of the hierarchy between Mx and the weak scale. We therefore usually 
allow superparticles as well as ordinary particles to be produced in X decays, as described by 
the minimal supersymmetric extension of the Standard Model (MSSM). 

Among the top-down scenarios, we have to distinguish between (at least) two classes of X 
particles: GUT particles (gauge or Higgs bosons of masses close to the GUT scale ~ 10^® GeV) 
could be trapped in topological defects (cosmic strings or monopoles) formed during a phase 
transition in the early Universe, and might be released when these defects annihilate or collapse 
m Another possibility is to consider metastable X particles HU that are distributed freely 
in the Universe; they might have been created at the end of inflation HH If their lifetime 
Tx exceeds the age of the Universe by about a factor 10^° ■ (10^^ GeV)/Mx, the X density 
required to explain the observed flux of UHEGR events also makes X a good cold Dark Matter 
candidate [nmn]. Depending on the model we are considering, X decay modes can be very 
different. For example, GUT gauge bosons will undergo two-body decays into all light (MS)SM 
particles with comparable branching ratios. On the other hand, if X decays only proceed 
through higher dimensional operators, X will usually undergo many-body decays. In another 
recently proposed “brane world” scenario ^3], Tx is exponentially suppressed by the relatively 
large distance between two branes, and X decays predominantly into light Higgs bosons and 
their superpartners. In any case, it is reasonable to assume that X will decay into some light 
particles contained in the spectrum of the MSSM. These primary decay products initiate parton 
cascades, the development of which can be predicted from the “known” interactions contained 
in the MSSM. 

Here we describe in detail the calculation of the spectrum of stable particles (protons, 
electrons, photons, three kinds of neutrinos, and lightest superparticles LSP) produced in X 
decays; a brief summary of these results has appeared in HU- The hrst calculations of this 
kind [S] used simple scaling fragmentation functions to describe the transition from partons to 
hadrons. Later analyses [nmn] used Monte Garlo programs to describe the cascade. However, 
since we can only expect to see a single particle from any given cascade, we only need to know 
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the one-particle inclusive decay spectrum of X. This is encoded in fragmentation functions; the 
evolution of the cascade corresponds to the scale dependence of these FFs, which is described 
by generalized DGLAP equations nn. In contrast to earlier analyses using this technique, 
which only include (SUSY) QCD [THl ITTH OlH 1^ . or at best a partial treatment of electroweak 
interactions OH] , we consider all gauge interactions as well as third generation Yukawa 
interactions; note that at energies above 10^° eV, due to the running of the coupling constants, 
all gauge interactions are of comparable strength. We will show that the inclusion of electroweak 
gauge interactions in the shower gives rise to a signihcant flux of very energetic photons and 
leptons, beyond the highest proton energies. Moreover, we carefully model decays of all unstable 
particles. As a result, we are for the hrst time able to fully account for the energy released in 
X decay. We cover all possible primary X decay modes, i.e. our results should be applicable 
to all models where physics at energies below Mx is described by the MSSM. 

Before going into technical details, we briefly outline the physics involved in the decay of 
an ultra-heavy X particle; it is summarized in £g. ^ By assumption its primary decay is into 
2 or more particles of the MSSM. These primary decay products will generally not be on-shell; 
instead, they have very large (time-like) virtualities, of order Mx- Each particle produced in 
the primary decay will therefore initiate a parton shower. The basic mechanism driving the 
shower development is the splitting of a virtual particle into two other particles with (much) 
smaller virtualities; the dynamics of this process is described by a set of splitting functions 
(SFs). As long as the virtuality is larger than the typical sparticle mass scale Msusy, all 
MSSM particles participate in this shower. At virtuality Msusy ~ 1 TeV the breaking of both 
supersymmetry and of SU{2) x U{1)y gauge invariance becomes important. All the massive 
superparticles that have been produced at this stage can now be considered to be on-shell, 
and will decay into Standard Model (SM) particles and the only (possibly) stable sparticle, 
the LSP. The same is true for the heavy SM particles, i.e. the top quarks and the massive 
bosons. However, the lighter quarks and gluons will continue a perturbative parton shower 
until they have reached either their on-shell mass scale or the typical scale of hadronization 
Qhad ~ 1 GeV. At this stage, strong interactions become non-perturbative, forcing partons to 
hadronize into colorless mesons or baryons. Finally, the unstable hadrons and leptons will also 
decay, and only the stable particles will remain. The spectra of these particles constitute the 
result of our calculation, which gives the UHEGR spectrum at the location of X decay. Of 
course, the spectrum on Earth might be modihed considerably due to propagation through the 
(extra)galactic medium [Hj; we will not address this issue here. 

Technically the shower development is described through fragmentation functions (FFs). 
The dependence of these functions on the virtuality is governed by the DGLAP evolution 
equations na extended to include the complete spectrum of the MSSM. All splitting functions 
needed in this calculation are collected in Appendix A. We numerically solved the evolution 
equations for the FFs of any particle of the MSSM into any other. At scale Msusy we applied 
unitary transformations to the FFs of the unbroken fields (“current eigenstates”) in order to 
obtain those of the physical particles (“mass eigenstates”); details are given in Appendix B. We 
then model the decays of all particles and superparticles with mass ~ Msusy, using the public 
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9 ^. 




(SUSY (hadronization) 
+ SU{2)®U(l) 
breaking) 


Figure 1: Schematic MSSM cascade for an initial squark with a virtuality Q ~ Mx- The 
full circles indicate decays of massive particles, in distinction to fragmentation vertices. The 
two vertical dashed lines separate different epochs of the evolution of the cascade: at virtuality 
Q > Msusy, all MSSM particles can be produced in fragmentation processes. Particles with 
mass of order Msusy decay at the hrst vertical line. For Msusy > Q > Qhad light QCD degrees 
of freedom still contribute to the perturbative evolution of the cascade. At the second vertical 
line, all partons hadronize, and unstable hadrons and leptons decay. See the text for further 
details. 


code ISASUSY j21j to compute the branching ratios for all allowed decays, for a given set of 
SUSY parameters. If R-parity is conserved, we obtain the hnal spectrum of the stable LSP 
at this step; the rest of the available energy is distributed between the SM particles. After a 
second perturbative cascade down to virtuality ~ max(mq, Qhad), the quarks and gluons will 
hadronize, as stated before. This non-perturbative phenomenon is parameterized in terms of 
“input” FFs. We use the results of ref.j22j, which are based on fits to LEP data. We paid special 
attention to the conservation of energy; this was not possible in previous studies, because of 
the incomplete treatment of the decays of particles with mass of order Msusy- We are able to 
check energy conservation at each step of the calculation, up to a numerical accuracy of a few 
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per mille. 

The remainder of this article is organized as follows. In sec. 2 we describe the technical 
aspects of the calculation. The derivation and solution of the evolution equations is outlined. 
We also check that our hnal results are not sensitive to the necessary extrapolation of the 
input FFs. Numerical results are presented in sec. 3. We give the energy fractions carried by 
the seven stable particles for any primary X decay product, and study the dependence of our 
results on the SUSY parameters. We finally describe our implementation of color coherence 
effects at small x using the modified leading log approximation (MLLA). Sec. 4 is devoted to 
a brief summary and conclusions. Technical details are delegated to a series of Appendices, 
giving the complete list of splitting functions, the unitary transformations from the interaction 
states to the physical states, our treatment of 2- and 3-body decays, parameterizations of the 
input FFs, and hnally a complete set of FFs obtained with our program for a given set of 
SUSY parameters (corresponding to a gaugino-like LSP with a low value of tan/? ~ 3.6 and 
AfsusY ~ 500 GeV). 


2 Technical aspects of the calculation 


In this section we describe how to calculate the spectra of stable particles produced in X decays: 
protons, electrons, photons, the three types of neutrinos and LSPs, and their antiparticles. Note 
that at most one out of the many particles produced in a typical X decay will be observed on 
Earth. This means that we cannot possibly measure any correlation between different particles 
in the shower; the energy spectra of the final stable particles are indeed the only measurable 
quantities. These spectra are given by the differential decay rates dTx/dEp, where P labels 
the stable particle we are interested in. This is a well-known problem in QCD, where parton 
showers were hrst studied. The resulting spectrum can be written in the form 


dxp 


E 


dT{X I) 
dxj 



( 1 ) 


where / labels the MSSM particles into which X can decay, and we have introduced the scaled 
energy variable x = 2E/Mx- dT{X —> I)/dxi depends on the phase space in a particular decay 
mode; for a two-body decay, dV{X —^ I)/dx[ oc 5(1 — xj). The convolution is dehned as 

f{z)®g{x/z)= f f{z)g(-) —. (2) 

Jx \Z J Z 

All the nontrivial physics is now contained in the fragmentation functions (FFs) Df {z,Q‘^). 
They encode the probability for a particle P to originate from the shower initiated by another 
particle I, where the latter has been produced with initial virtuality Q. This implies the 
“boundary conditions” 

= ^/■ ^(1 - (3) 

which simply say that an on-shell particle cannot participate in the shower any more. As 
already explained in the Introduction, for Q > Msusy all MSSM particles J are active in the 
shower, and thus have to be included in the list of “fragmentation products”. 
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The evolution of the FFs with increasing virtuality is described by the well-known DGLAP 
equations HU- In the next two subsections we discuss these evolution equations, and their 
solution, in more detail. We first only include strong (SUSY-QCD) interactions. However, at 
energies above 10^° eV all gauge interactions are of comparable strength. The same is true for 
interactions due to the Yukawa coupling of the top quark, and possibly also for those of the 
bottom quark and tan lepton. In a second step we therefore extend the evolution equations to 
include these six different interactions. We then describe the decays of heavy (s)particles, which 
happen at virtuality Q = Msvsy ~ 1 TeV. At Q < Msusy only QCD interactions need to be 
included, greatly simplifying the treatment of the evolution equations in this domain. Finally, 
we describe the nonperturbative hadronization, and the weak decays of unstable hadrons and 
leptons. 


2.1 Evolution equations in QCD and SUSY-QCD 

For convenience, we review here the DGLAP evolution equations in ordinary QGD. As already 
noted, the FF Dp{x,Q‘^) of a parton (quark or gluon) p into a particle (parton or hadron) P 
describes the probability of fragmentation of p into P carrying energy Ep = xEp at a virtuality 
scale Q. If P is itself a parton, the FF has to obey the boundary condition ( 0 . However, if 
F is a hadron, the x—dependence of the FF cannot be computed perturbatively; it is usually 
derived from fits to experimental data. Perturbation theory does predict the dependence of the 
FFs on the virtuality Q: it is described by a set of coupled integro-differential equations. In 
leading order (LO), these QGD DGLAP evolution equations can be written as |2ni : 


dD!;,{x,Q^) 
d log((32) 
dD;{x,Q^) 
d log((32) 



where 05 is the running QGD coupling constant, E is the number of active flavors (i.e. the 
number of Dirac quarks whose mass is lower than Q), and i labels the quarks and antiquarks.1 
The convolution has been defined in eq. ©• The physical content of these equations can be 
understood as follows. A virtual quark g* can reduce its virtuality by emitting a gluon; the 
hnal state then contains a quark and a gluon. Either of these partons (with reduced virtuality) 
can fragment into the desired particle F; this explains the occurrence of two terms in the hrst 
eq. ©• Analogously, a gluon can either split into two gluons, or into a quark-antiquark pair, 
giving rise to the two terms in the second eq.(j 3 ]). 

These partonic branching processes are described by the splitting functions (SFs) Pp^p^^x), 
for parton Pi splitting into parton p 2 , where x = Fp^/Fp^. As already noted, in pure QGD there 

tNote that the DGLAP equations given here are the time-like ones, which describe the evolution of fragmen¬ 
tation functions. In leading order they differ from the space-like DGLAP equation (describing the evolution 
of distribution functions of partons inside hadrons) only through a transposition of the matrix of the splitting 
functions. 
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are only three such processes: gluon emission off a quark or gluon, and gluon splitting into a qq 
pair. The hrst of these processes gives rise to both SFs appearing in the hrst eq. dH) ; momentum 
conservation then implies Pqq{x) = PgqiX — x), for a: 7 ^ 1. Similarly, Pgg{x) = Pggi)- — x) and 
Pqg{x) = Pqg{^ — x) follows from the symmetry of the hnal states resulting from the splitting 
of a gluon. Special care must be taken as a; —> 1. Here one encounters infrared singularities, 
which cancel against virtual quantum corrections. The physical result of this cancellation is 
that the energy of the fragmenting parton p is conserved, which requires 

Y. [ xD^{x,Q^) = 1 \/p,Q^. (5) 

p ^ 0 


This can be ensured, if 


f dxxY.Py'vi^) =0 '^P- 


Note that these integrals must give zero (rather than one), since eqs. 
of the FFs. The explicit form of the QCD SFs is [T7j : 


( 6 ) 


only describe the change 


Pggix) 


4 /1 + a:^ 
3 


X 


Pgqix) 


4 1 + (1 — x)'^ 

3 X ’ 


2 


Pgg{x) = 6 


X 


X 


+ a;(l — x) + 


X 


(1 -x). 


+ i(i-x)(F-n 


F 

18 


(7) 


The “+” distribution, which results from the cancellation of x - 
above, is dehned as: 

/ f{x)g{x)+dx = f [/(x) - f{l)]g{x) 

JO JO 


1 divergences as outlined 


( 8 ) 


while 5 '(x)+ = g{x) for x 7 ^ 1. Finally, the scale dependence of is described by the 

following solution of the relevant renormalization group equation (RGE): 


asm = 


2tiB 

iTf 


(9) 


where A ~ 200 MeV is the QCD scale parameter, and B = 6/(33 — 2F). 

Note that eqs. © list different FFs for all (anti)quark flavors g*. At hrst sight it thus seems 
that one has to deal with a system of 2F + 1 coupled equations. In practice the situation can 
be simplihed considerably by using the linearity of the evolution equations. This implies 


Dpx,Q^) = Y.Dt(^,Q'\Ql)»Dp-,Ql), 

P' ^ 


( 10 ) 
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where the generalized FFs again obey the evolntion eqnations iD- Moreover, they satisfy 
the boundary conditions D^' {x, Ql, Ql) = 5(1 — x)6^' at some convenient value of Qq < Q. 
The D thus describe the purely perturbative evolution of the shower between virtualities Q 
and Qo- This ansatz simplihes our task, since all quark flavors have exactly the same strong 
interactions, i.e. we can use the same D^. for all quarks g* with < Qq. Moreover, we only 
have to distinguish three different cases for p' : qi, qj with j ^ i, and g. All flavor dependence 
is then described by the Dp,(x,Qq); for sufficiently small Qq, these can be taken directly from 
hts to experimental data. If we make the additional simplifying assumption that all quarks 
and antiquarks are produced with equal probability in primary X decays, we effectively only 
have to introduce two generalized FFs D for a given particle P, one for the fragmentation of 
gluons and one for the fragmentation of any quark. In other words, in pure QCD we only need 
to solve a system of two coupled equations. 

The introduction of squarks g* and gluinos g, i.e. the extension to SUSY-QCD, requires the 
introduction of FFs D^., D?. This gives rise to new SFs, describing the emission of a gluon 
by a squark or gluino, as well as splittings of the type g* —>• g^g, g* —>• g^g and g —>• g^gj. We 
thus see that any of the four types of partons (g*, g*, g, g) can split into any (other) parton. The 
complete set of evolution equations thus contains 16 SFs |2Z1, which we collect in Appendix A. 
The presence of new particles with SU{3) interactions also modihes the running of as- One 
can still use eq.©, but now Bsusy = 2/(9 — F). 

2.2 Evolution equations in the MSSM 

We now extend our discussion of the evolution equations to the full MSSM. We already saw in 
the Introduction that superparticles can only be active in the shower evolution at virtualities 
Q > Msusy 1 TeV. This means that the supersymmetric part of the shower evolution can 
be described in terms of generalized FFs Dj satisfying the boundary condition 

i)/(a;, M|usy) ^susy) = ~ ^)i (H) 

where I and J label any (s)particle contained in the MSSM. Note that ea. (ITT|l differs from eq.Q 
since the former is valid for all particles in the MSSM, including light partons. According to 
the discussion following eq. ca we only have to consider those particles to be distinct that have 
different interactions. We include all gauge interactions in this part of the shower evolution, 
as well as the Yukawa interactions of third generation (s)fermions and Higgs bosons, but we 
ignore hrst and second generation Yukawa couplings, as well as all interactions between different 
generations. This immediately implies that we do not need to distinguish between hrst and 
second generation particles. Moreover, we ignore CP violation, which means that we need not 
distinguish between particles and antiparticles. Finally, the electroweak SU{2) symmetry can 
be taken to be exact at virtuality Q > Msusy, he. we need not distinguish between members 
of the same SU{2) multiplet. This is analogous to ordinary QCD, where one does not need 
to introduce different FFs for quarks with different colors. Our assumption implies that X is 
an SU{2) singlet. Had we allowed |221 ^ to transform nontrivially under SU{2), the SU{2) 


splitting functions would have to be modified I2H|.* 

Altogether we therefore need to treat 30 distinct particles: six quarks 
four leptons e^j, r^, r/j, three gauge bosons B,W,g, two Higgs bosons Hi,H 2 , and all their 
superpartners; Hi couples to down-type quarks and leptons, while H 2 couples to up-type 
quarks. Note that a “particle” often really describes the contribution of several particles which 
are indistinguishable by our criteria. For example, the “quark” un stands for all charge—2/3 
right-handed quarks and antiquarks of the two first generations, i.e. ur, cr and their antipar¬ 
ticles ur^ Cr. This can be expressed formally as = {Dur + ^cr + -^C?) /^’ where in 

our approximation the four terms in the sum are all identical to each other after the final state 
P has been summed over particle and antiparticle.^ Similarly, qr stands as initial particle for 
an average over the two SU{2) quark doublets of the two first generations (m^, d^) and {cr, sr), 
and their antiparticles. Note that all group indices of the particle in question are summed over. 
In the usual case of QCD this only includes summation over color indices, but in our case it 
includes summation over SU{2) indices, since SU{2) is (effectively) conserved at energies above 
TfsUSY- 

Let us hrst discuss the scale dependence of the six coupling constants that can affect the 
shower evolution significantly at scales Q > Msusy- These are the three gauge couplings gy, 
g 2 and gs, which are related to the corresponding “fine structure constants” through = 
gf/{4n), i G {1/2,5'}. Moreover, the third generation Yukawa couplings are proportional to 
the masses of third generation quarks or leptons: 


yt 

Vh 

Vr 


gmt 

\/2 mw sin (3 ’ 
9 mb 

\/2 my/ cos (3 ’ 
gmr 

■\/2 mw cos (3 ’ 


( 12 ) 


where tan/3 = {H 2 )/{Hi). The couplings i/b and yr are only significant if tan/3 S> 1. Note 
that in many models, values tan/3 ~ mtimt)/mbimt) — 60 are possible, in which case yb and 

tA minor caveat might be in order here. We assume the electroweak part of the evolution of the FFs to be 
described by the equations for an unbroken SU (2) x U{1)y theory. It seems physically reasonable that this should 
be a good approximation as long as both the energy and the virtuality of the gauge bosons is large compared to 
their masses, which is the case in our calculation. Indeed, existing loop calculations m confirm that the leading 
logarithmic corrections, which we treat here, are the same in the broken and unbroken theory, as long as all scales 
are large compared to Mw- However, these calculations all deal with the electroweak equivalent of structure 
functions, rather than fragmentation functions. Strictly speaking our assumption, although reasonable, has 
therefore not yet been backed up by an explicit calculation that includes massive gauge bosons. 

consistent interpretation of, e.g., ur as a “particle” requires that ur stands for the average of ur^ cr etc. 
when Ur appears as lower index of a generalized FF, as described in the text. However, ur stands for the sum 
of Ur, Cr etc. when ur is an upper index of a D. With this definition, we have .D“«(x, Mgugy) = <5(1 —x). This 
interpretation also fixes certain multiplicity factors in the DGLAP equations, as detailed in Appendix A. This 
treatment is only possible if X has equal branching ratio into ur, cr etc. However, we expect the differences 
between decays into first or second generation quarks to be very small even in models where these branching 
ratios are not the same. 
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Ur are comparable in magnitude to gs and ( 72 , respectively. The LO RGBs for these six MSSM 
couplings are |BD|: 


dt 

dg2 

dt 

dgs 

dt 

dyt 

dt 

dVh 

dt 

dVr 

dt 


11 ^ 

167r2 ’ 

92 

o 9s 




yt 

167r2 

yt 

yr 

167r2 


(^yt + yl - -jdY - 3^2 - y, 
+ yt +yl- ^ 9 y - 3^2 - Yds 
{^yt “ 1 “ ~ ^9y ~ ^9 


(13) 


where t = log ^ parameterizes the logarithm of the virtuality, and Qo is an arbitrary scale 
where the numerical values of these couplings constants are “known” (in case of the Yukawa 
couplings, up to the dependence on tan/?). As well known [211, given their values measured at 
Qo — 100 GeV eas. dTTH) predict the three gauge couplings to unify at scale Mqut — 2-10^® GeV, 
i.e. S'KMqut) = 5 ' 2 (^gut) = 59y{^gvt)/^ — 0.52, where the Glebsch-Gordon factor of 5/3 is 
predicted by most simple unihed groups, e.g. SU{5) or S'O(IO). We solved these equations by 
the Runge-Kutta method; of course, the RGBs for the gauge couplings can trivially be solved 
analytically, but the additional numerical effort required by including eqs. m in the set of 
coupled differential equations that need to be solved numerically is negligible. 

The main numerical effort lies in the solution of the system of 30 coupled DGLAP equations, 
which are of the form: 


dbj 

d\og{Qb 


(a:, Q^, Mlugy) — y/ 

K 


gki{Q^ 


-Pki{z) 0 Il^( —, Q , Afgygy) , 

Z 


(14) 


where I,J,K run over all the 30 particles, and axiiQb = 91:11 is the (running) coupling 
constant associated with the corresponding vertex; note that at this stage we are using interac¬ 
tion (or current) eigenstates to describe the spectrum. Generically denoting particles with spin 
1, 1/2 and 0 as V, F and S (for vector, fermion and scalar), we have to consider* branching 
processes of the kind V ^ VV, V ^ FF, V ^ SS, F ^ FV, F ^ FS, S ^ SV and 
S —>■ FF. All these branching processes already occur in SUSY-QGD. The splitting functions 
can thus essentially be read off from the results of ref. ^7], after correcting for different group 
[color and/or SU{2)] and multiplicity factors. The coefficients of the 5(1 — x) terms in diagonal 
SBs can be hxed using the momentum conservation constraint in the form ©; note that these 

*We do not need to consider S SS, since the corresponding dimensionful coupling is 0{Ms\jsy) <C Q in 
this domain, i.e. these processes are much slower than the relevant time scale l/Q. 
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constraints have to be satisfied for each of the six interactions separately. The explicit form of 
the complete set of MSSM SFs Pki{x) is given in Appendix A. 

We solved these eqnations nnmerically using the Runge-Kutta method. To that end the 
FFs were represented as cubic splines, using 50 points which were distributed equally on a 
logarithmic scale in x for 10“^ ^ x < 0.5, and 50 additional points distributed equally in 
log(l — x) for 0.5 < X < 1 — 10“^. Starting from the boundary conditions^ (j1 1 jl . we arrive at 
the 30 X 30 generalized fragmentation functions at virtuality Q = Mx- Here we assume that 
the evolution equations describe the perturbative cascade at these energies correctly. We will 
comment on the limitations of our treatment at the end of this Section. 

2.3 Evolution of the cascade below Q = 1 TeV 

Here we would like to describe the physics at scales at and below Msvsy- the breaking of both 
supersymmetry and SU{2) ® U{1) symmetry, the decay of unstable (s)particles with masses of 
order Msusy, the pure QCD shower evolution down to Qhad, the non-perturbative hadroniza- 
tion of quarks and gluons, and finally the weak decays of unstable leptons and hadrons. For 
simplicity we assume that all superparticles, the top quark as well as the W, Z and Higgs bosons 
all decouple from the shower and decay at the same scale Msvsy — 1 TeV. The fragmentation 
of b and c quarks is treated using the boundary condition © at their respective mass scales 
of 5 and 1.5 GeV, while the nonperturbative hadronization of all other partons takes place at 
Qhad = 1 GeV. 

At Q = Msusy we break both Supersymmetry and SU{2) ® U{1). All (s)particles acquire 
their masses in this process, and in many cases mix to give the mass eigenstates. This means 
that we have to switch from a description of the particle spectrum in terms of current eigen¬ 
states to a description in terms of physical mass eigenstates. This is accomplished by unitary 
transformations of the type* * 

Cf = E . (15) 

J 

Unitarity requires Es |csjp = Ej |csjp = 1, if the current state J has the same number of 
degrees of freedom as the physical state S. This is often not the case in the usual convention; 
then some care has to be taken in writing down the |c 5 jp, see Appendix B. We use the 
following physical particles: u ,d ,b, s ,c ,t quarks and e , ,t leptons now have both left- and 
right-handed components, i.e. they have twice as many degrees of freedom as the corresponding 
states with fixed chirality. The neutrinos remain unchanged, since we ignore the interactions of 
right-handed neutrinos. The gluons also remain unchanged, since SU{3) remains exact below 
Msvsy- The electroweak gauge sector of the SM is described by W := + W~, Z and 

7; note that the massive gauge bosons absorb the Goldstone modes of the Higgs sector, and 
hence receive corresponding contributions in eq. (USD. The Higgs sector consists of two charged 

^Technically, these 5—functions are represented by narrow Gaussians centered at a: = 1, normalized to give 
unity after integration over x < 1. 

*Note that the squares of the coefficients csj appear in eq. m, since the FFs describe probabilities, which 
are related to the square of the wave functions of the particles in question. 
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Higgs bosons (described by if = + H~) and the three neutral ones and 

the neutral Higgs bosons are described by real fields, which contain a single degree of freedom. 
In the SUSY part of the spectrum, the gluino g as well as the first and second generation 
sfermions ul,r, cl,r and sneutrinos remain unchanged (but ur and (Ir, etc., are now 

distinguishable). The SU{2) singlets and doublets of third generation charged sfermions mix 
to form mass eigenstates U, ^ 2 , b 2 , fi, f 2 - Similarly, the two Dirac charginos xf ^md yfare 

mixtures of charged higgsinos and winos, and the four Majorana neutralinos y®; Xs; X 4 ) 
order of increasing masses, are mixtures of neutral higgsinos, winos and binos. 

The numerical values of many of the crj depend on the parameters describing the breaking 
of supersymmetry. We choose four different sets of parameters, which describe typical regions 
of the parameter space, in order to study the impact of the details of SUSY breaking on the final 
spectra. We take two fairly extreme values of tan(/9) = 3.6 and 48, and two sets of dimensionful 
parameters corresponding to higgsino-like and gaugino-like states yf, Xi and X 2 - We used the 
software IS ASUS Y [21] to compute the mass spectrum and the mixing angles of the sparticles 
and Higgses for a given set of SUSY parameters. 

Having computed the spectrum of physical (massive) particles, we have to treat the decay 
of all unstable particles with mass near Msusy- Since we assumed i?—parity to be conserved, 
the lightest supersymmetric particle (LSP) is stable. In our four scenarios (as in most of 
parameter space) the LSP is the lightest neutralino Xi- The end products of these decays 
are thus light SM particles and LSPs. Note that decays of heavy sparticles often proceed 
via a cascade, where the LSP is produced only in the second, third or even fourth step, e.g. 
g —> uul —> udxt ude'^h'eXi- order to model these decays we again use ISASUSY, 
which computes the branching ratios for all allowed tree-level 2- and 3-body decay modes of 
the unstable sparticles, of the top quark and of the Higgs bosons. Together with the known 
branching ratios of the W and Z bosons, this allows us to compute the spectra of the SM 
particles and the LSP after all decays, by convoluting the spectra of the decaying particles 
with the energy distributions calculated for 2- or 3-body decays. The total generalized FF of 
any MSSM current eigenstate I into a light or stable physical particle s (quark, gluon, lepton, 
photon or LSP) is then 

D^j = Df=^ + Y.Df®PsS, (16) 

where Pss describes the spectrum of s in the decay S' —>• s. We compute these spectra from 
phase space, including all mass effects, but we didn’t include the matrix elements. The spectra 
for each decay mode of the heavy particle S are normalized to give the correct branching ratio, 
as computed by ISA JET. As far as LSPs are concerned, ea. ([TT)|l already gives the final result, 
i.e. If s is a lepton or photon, eq. (USD describes the FF at all virtualities between 

Msvsy and nib = 5 GeV. 

As we will see shortly, in some cases two-body decays can lead to sharp edges in the FFs at 
intermediate values of x. This can happen if the primary decay product is a massive particle 
with only weak interactions. In that case a substantial fraction of the initial J—peak at a; = 1 
survives even after the evolution; convolution of this J—peak with a two-body decay distribution 
leads to a fiat x distribution of the decay products between some Xmin and a^max- An accurate 
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description of these contributions to the FFs sometimes requires the introduction of additional 
points near Xmin and/or Xmax in the splines describing these FFs. 

The perturbative evolution in the QCD sector does not stop at Msusy, but continues until 
virtuality Qq = max(mq, Qhad)- This part can be treated by introducing generalized FFs as 
in eq. (uni), where (p, p') G {u, d, s, c, b, gj are light QCD partons. We use once more the DGLAP 
evolution equations, but this time for pure QCD, evolving these generalized FFs between Qo 
and Msusy- The generalized partonic FFs between Qo and Mx can then be computed through 
one more convolution: 


D^{x, Mx, Qo) — ^ Dj {z, Mx, Mgjjgy) ® Dp,( , Mgjjgy, Qo) • 


(17) 


The total partonic FFs at Mx can finally be computed through eq. m by using known 
“input FFs”. They describe the non-perturbative hadronization of quarks and gluons into 
mesons and baryons, which happens at Q = Qo. These FFs D/(x, Qq), where i G {u, d, s, c, b, g} 
and h represents a hadron, can be obtained directly from a £t to (e.g.) LEP data. We 
used the results of where the FFs of a quark or gluon into protons, neutrons, pions and 
kaons (or more exactly the sum over particles and antiparticles) are parameterized in the form 
iVx“(l — x)^. 

The original form [23j of these functions is only valid down to x = 0.1. Kinematic and color 
coherence effects, which are not included in the usual DGLAP framework, become important 
[B2]) at X ^ (Q/Qhad) ~ 0.1, where in the second step we have used the LEP energy scale 
Q ~ 100 GeV. For Q ~ Mx ~ 10^® GeV these effects become large only for x ^ 10“®; they can 
thus safely be ignored for many (but not all; see below) applications. In we therefore chose 
a rather simple extrapolation of the functions given in I2SI towards small x. Our default choice 
was a Nx~°‘ parameterization; N and a' were computed by requiring the continuity of this 
parameterization with the FFs of |2Sj at some Xq — 0.1, energy conservation, and, as additional 
constraint, an identical power law behavior at small x (i.e. identical a') for all the FFs of a 
given quark into the different hadrons. This last assumption was motivated by the fact that we 
obtain such an identical power law at small x during the perturbative part of the cascade, and 
by the well accepted LPHD hypothesis (Local Parton-Hadron Duality) jSSI, which postulates 
a local proportionality in phase space between the spectra of partons and hadrons. We chose 
different xq for each initial parton in such a way that we obtain a' between 0 and 2; the upper 
bound on a' follows from energy conservation (the energy integral Jq dxxD(x) has to be finite). 

In order to check the consistency of this parameterization, we used another functional form 
with three free parameters: D{x) = ax““ + 61ogx + c, a > 0. This allowed us to freely choose 
a', keeping the same assumptions about continuity etc. as above. This enabled us to compare 
two extreme values of a', namely 0.5 and 1.4. The hrst is the smallest value compatible 
with a > 0, while the second approximates the small—x behavior of the perturbative QGD 
evolution between 1 GeV and 1 TeV; requiring a' < 1.4 thus ensures that this perturbative 
evolution dominates the behavior of the FFs at small x. Note also that the perhaps most 
plausible value, a' ~ 1 (which corresponds to a flat distribution of particles in rapidity when 
perturbative effects are ignored) is comfortably bracketed by these limiting values. In fig. |21we 
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Figure 2: Effect of varying the low—a; extrapolation of the input FFs on the hnal FFs 
and D'^ . See the text for further explanations. 


plot the hnal result at small x for different FFs with these two extreme parameterizations, after 
convolution with the perturbative FFs. As can be seen, the effect of varying a' is very small 
once energy conservation is imposed. This indicates that our hnal results are not sensitive to 
the necessary small—x extrapolation of the input FFs.^ The main uncertainty at moderately 
small X (10“^^x^0.1) will then come from perturbative higher order corrections, which might 
be quite signihcant in this range. 

Unfortunately, we were not able to perform a complete NLO analysis, for the following 
reasons. Beyond leading order the SFs for space-like and time-like processes are no longer 
identical Isa. Already at next-to-leading order (NLO) the time-like SFs have a rather bad 
behavior at small x, with a negative leading term — in Pqq. This term is tempered in the 
hnal spectra (which have to be positive) by the convolution occuring in the DGLAP equations, 
as well as by the convolution of the FFs with NLO “coefficient functions” which modify the 
basic relation o once higher order corrections are included. Note that the FFs, SFs and NLO 
coefficient functions are scheme dependent; worse, the coefficient functions are also process- 
dependent, i.e. they will depend on the spins of X and its primary decay products. NLO 
results are known for the classical processes occuring in pure (non-supersymmetric) QCD, but 
they are not available for most of the processes we are interested in. Moreover, in cases where 
they are known, these coefficient functions often contain the most important part of the NLO 
correction, rendering useless any attempt to give a partial result by only including NLO terms 
in the SFs. We conclude that it might be possible and interesting to carry out a full NLO 

’^However, the original FFs of ref.[21] should not be used on the whole range [10“”^, 1], since they violate 
energy conservation badly, leading to over-production of particles at small x. 
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analysis in the pure QCD case, but this is not possible in the more interesting supersymmetric 
case using available results. Note that part of the perturbative NLO effects are absorbed in the 
input FFs, through their fit to experimental data. At very small x, NLO effects just give the 
leading “color coherence” corrections, which are resummed analytically in the MLLA formula, 
as will be discussed in Sec. 3.4. 

Finally, having computed the spectra of long-lived hadrons and leptons, we still need to 
treat weak decays of unstable particles, in order to obtain the final spectra of protons, electrons, 
photons and the three types of neutrinos. This is again done using the formalism of eq. (fT^ . We 
limit ourselves to 2- and 3-body decays, considering the 4-body decays of the r to be cascades 
of 2-body decays. As before, we compute the decay functions Psh for FT — > s decays from phase 
space only, and we ignore decays with branching ratio smaller than 1%. We then renormalize the 
branching ratios of the decays we do include, so that we maintain energy conservation. We also 
explicitly treated the leptonic part of the semi-leptonic decays of b— and c—flavored hadrons, 
which are evidently not included in the FFs of |2Sj. We used the Peterson parameterization for 
non-perturbative heavy quark fragmentation |2Sj, and then treated the semi-leptonic decays 
in the spectator model (i.e. using the same spectra as for free quark decays, with rric = 1.5 
GeV and mb = 4.5 GeV). Details of our treatment of decays are given in Appendix G. 

3 Results and analysis 

3.1 General features of the final fiuxes 

Some FFs computed with our code have already been presented in ref. ng. In particular, 
we showed that including supersymmetric particles in the QGD shower evolution signihcantly 
softens the FFs of strongly interacting (s)particles. Usually FFs are multiplied with to 
allow a more direct comparison with the experimental UHEGR spectra, which are multiplied 
with to make them (approximately) flat. In this normalization, including supersymmetry 
reduces the height of the peak of the FFs of initial quarks and gluons by approximately a 
factor of two. Including electroweak interactions reduces this peak by another 20% or so. More 
importantly, it leads to a very energetic component of the final FFs into leptons and especially 
photons, if the primary particle has electroweak interactions. This is due to the emission of 
very energetic electroweak gauge bosons early in the parton shower; some of them will decay 
or fragment (partly into leptonic final states), while others survive as photons after the unitary 
transformation m- In this article we focus on a detailed description of X decays computed 
according to the “state of the art”, i.e. with all interactions included. 

A fairly complete set of results of our code for a given set of SUSY parameters is given in 
Appendix E. Here we assumed similar masses for all sfermions, higgsinos, heavy Higgs bosons 
and gluinos, ~ ~ ~ /x ~ 500 GeV; this leads to a gaugino-like LSP, since we assume 

“gaugino mass unification”, i.e. 6m^ ~ 3m^y ~ mg. We also choose a small value for the ratio 
of vevs, tan/3 = 3.6. We see that the final spectra depend sensitively on the primary X decay 
products especially in the large x region. This strong dependence on the unknown primary 
X decay mode(s) should be kept in mind when one is trying to quantitatively test “top-down” 
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models. Nevertheless, we can make a few general statements abont these resnlts. To that end 
we hrst analyze ratios of FFs of the different stable particles divided by the FF of the same 
initial particle into protons. Recall that these FFs directly represent the flnx at sonrce if X 
nndergoes two-body decay. 

Taking the ratios of the different FFs renders some featnres more evident, as can be seen 
from figs. El and m First of all, in the low x region most FFs show the same power law behavior, 
and the ratios become qnite independent of the initial particle. The exceptions are the FFs 
into the LSP and Ur- This comes from the fact that the LSP flnx as well as most of the 
Ur are prodnced in the pertnrbative cascade above 1 TeV and in the following decays of the 
heavy particles of the spectrnm; they receive no contribntion from the decays of light hadrons, 
althongh the Ur flnx receives a minor contribntion from the decay of 6 —flavored hadrons. In 
contrast, at low x the flnxes of Ue, e and 7 all dominantly originate from the decays of 
light hadrons, in particnlar of charged or nentral pions; we saw in £g. |21 that the shape of the 
light hadron spectrnm at small x is essentially determined by the pertnrbative QCD evolntion, 
i.e. is independent of the initial particle I. In the region x ^ 0.01 we thns predict FFs into 
and 7 to be approximately 3 to 4 times larger than the FF into protons, while the FFs into 
electrons and Ug are aronnd twice the FF into protons. The FFs into LSP and Ur are five to 
20 times smaller than the one into protons. Note that the LSP flnx at small x from an initial 
particle is almost the same as that from its snperpartner. It is determined completely by the 
MSSM cascade, i.e. by the snpersymmetric DGLAP eqnations, and is almost independent of 
details of the snpersymmetric spectrnm. However, even at x = 0.01 the FF into the LSP does 
retain some sensitivity to the start of the cascade, i.e. to the initial particle I and hence to the 
primary X decay mode(s). 

At larger valnes of x the ratios of the FFs depend more and more strongly on the initial 
particle. As a; —>■ 1 the proton flnx is always orders of magnitnde smaller than the flnxes of all 
other stable particles. One reason is that the proton is a composite particle, i.e. its FF contains 
a convolntion with a non-pertnrbative factor which falls as a power of 1 — a; at large x. Even 
before this convolntion the flnx of partons (qnarks and glnons) that can give rise to protons is 
snppressed at large x dne to copions emission of (soft) glnons, whereas the FFs into leptons, 
photons and LSPs can remain large at large x. If the progenitor I of the cascade is a strongly 
interacting snperparticle, at large x the FF into the LSP always dominates over the other FFs. 
For an initial qnark or glnon, the flnx of 7 (which is the second after LSP for a sqnark or 
glnino) will dominate at large x. On the contrary, in the case of an initial lepton, PF, B or 
Hi, the strongest flnxes will be leptonic ones, the exact order depending of the initial particle. 
Moreover, for an initial (s)lepton, the flnxes will be significantly higher at high x (and hence 
smaller at low x, becanse of energy conservation) than for strongly interacting (snper)particles 
or Higgs bosons. Finally, an initial B or B has a 5—peak at a; = 1 (not visible in the fignres) in 
and respectively, in addition to a smooth component that vanishes as a; —> 1. This 

behavior reflects the inability of R or R to radiate a boson, i.e. there are no splitting processes 
R^R + AorR^R + X. 
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Figure 4: As in fig. El but for initial first or second generation SU{2) doublet lepton, I = 
(top) or slepton, I = II (bottom). 
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Table 1: Energy fractions Jq dxxD^{x, M\) carried by the stable particles p at the end of the 
cascade, for initial (s)quarks of the lst/2nd and 3rd generations. We took Mx = 10^® GeV and 
a sparticle spectrum with gaugino-like LSP, as described in Sec. 3.1. 


3.2 Energy distribution between the final stable particles 

In the following tables we show the total energy carried per each type of particle at the end of 
the cascade, depending on the progenitor of the cascade, for the same set of SUSY parameters 
as in Sec. 3.1. As stated earlier, we are able to verify energy conservation up to at most a few 
per mille at each step of the cascade, including its very end. We see that the “lost” energy is 
somewhat larger for (s)squarks, gluons and gluinos than for (s)leptons. This is due to numerical 
artefacts. The biggest numerical uncertainties arise from the Runge-Kutta method.^ 

Note that even for an initial quark or gluon, more than 35% of the energy is carried by 
the electromagnetic channels (electrons plus photons), while neutrinos carry about 40%; in 
this case most of these fluxes originate from the decays of light hadrons, chiefly pions. The 
corresponding numbers for superparticles are slightly smaller, the difference being made up by 
the increased energy fraction carried by the LSP (at large x)] an initial SU{2) singlet squark 
leads to a higher energy fraction in LSPs, since SU{2) singlet sfermions usually decay directly 
into the LSP, which is Bino-like for our choice of parameters, whereas SU{2) doublet sfermions 
preferentially decay via a cascade involving X2 xf- 

Lepton-induced showers have a far smaller photon component, but now an even larger 
fraction of the energy is carried by electrons and/or neutrinos, while protons carry at most 2% 

^■For practical reasons, we used a fixed virtuality step in this algorithm, which we had to keep reasonably 
large, the whole program being already quite time-consuming. In the worst cases, our choice of the virtuality 
step leads to errors of the order of a few per mille; such a precision is certainly sufficient for our purposes. 
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Table 2: Energy fractions carried by the stable particles at the end of the cascade, for initial 
(s)leptons of the lst/2nd and 3rd generations. Parameters are as in Table 1. 
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Table 3: Energy fractions carried by the stable particles at the end of the cascade, for initial 
bosons and bosinos. Parameters are as in Table 1. 
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of the primary’s energy. In this case the difference between an initial particle and its superpart¬ 
ner is much larger than in case of strongly interacting particles, since a much higher fraction of 
an initial slepton’s energy goes into LSPs, due to the reduced perturbative shower and shorter 
superparticle decay cascades. This also explains why more than 70% of the energy of an initial 
B [B) goes into photons (LSPs). The energy fractions for an initial SU{2) gauge or Higgs boson 
resemble those for a quark (with the exception of an increased z/,- component, which is however 
washed out by neutrino oscillations), although the shapes of the corresponding FFs differ quite 
dramatically. The energy fraction carried by protons is always quite small. Pions are created 
much more abundantly in the non-perturbative hadronization, and decay into leptons (2/3) 
and photons (1/3). As noted earlier, this explains the regularity and the features of the small 
X behavior. 

3.3 Dependence on SUSY parameters 

As stated in ra. the general features of our results described above depend very little on the 
set of SUSY parameters we are using. Here we give a more precise analysis of the influence 
of different parameters describing the SUSY spectrum. As usual we present our results as 
x^-D^{x, Mx). The multiplication with the third power of the energy leads to an approximately 
flat cosmic ray spectrum for E ^ 10^° GeV jHl- In our case it suppresses the small—a: region, 
leading to maxima in the curves at x between 0.1 and 1. 

We first studied the dependence of our results on the overall SUSY mass scale, by comparing 
results for two different IS ASUS Y input mass scales for scalars and gluinos: Msusy ~ 500 GeV 
and 1 TeV. As expected, this change has almost no impact on the final results, since the 
details of the decay chains of heavy (s)particles will depend mostly on the relative ordering 
of the (s)particle spectrum (e.g. allowing or preventing some decay modes), rather than on 
their absolute mass scale. Moreover, a factor 2 or 3 in the scale where the MSSM evolution is 
terminated does not change the FFs much, since the DGLAP equations describe an evolution 
which is only logarithmic in the virtuality. 

Next we compared two rather extreme values of tan/3, namely 3.6 and 48, leaving all dimen¬ 
sionful parameters at the weak scale unchanged. Once again the effect is very small. The only 
visible difference occurs for initial Hi and ifi, where the increase of tan/3 produces more Vr at 
large x, as can be seen in fig. El However, flavor oscillations will essentially average the three 
neutrino fluxes between source and detector, so we expect very little direct dependence of mea¬ 
surable quantities on tan /3. The main remaining effect is an increase of the overall multiplicity 
by 30% for an initial Hi or Hi in case of large tan /3, due to the increased shower activity 
from the much larger bottom Yukawa coupling. However, the situation could be different in 
more constrained models, where the spectrum is described by a few soft breaking parameters 
specified at some high energy scale. In this case a change of tan /3 generally changes the sparticle 
and Higgs spectrum, and can also greatly modify some branching ratios. 

In order to get a feeling for how the various FFs depend on the relative ordering of the 
dimensionful parameters describing the SUSY spectrum, we investigated two rather extreme 
cases. They resemble two qualitatively different regions of parameter space in the minimal su- 
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pergravity (mSUGRA or CMSSM) model where the thermal LSP relic density is acceptably 
small In the hrst scenario the LSP Xi hns small mass splitting to the lightest stan, 

Ti. We took the following valnes for the relevant soft breaking parameters: rriq ~ rrig = 1 
TeV for all sqnarks, = 250 GeV for all SU{2) donblet sleptons, mi^ ~ 200 GeV for 
I = e, fi bnt rednced rrif^ so that = m^o + 13 GeV = 163 GeV; note that in mSUGRA one 
needs large mass splitting between sqnarks and sleptons if the LSP mass is to be close to the 
Ti mass. The physical sfermion masses receive additional contribntions from SU{2) x U{1)y 
symmetry breaking, and, in case of the third generation, from mixing between singlet and 
donblet sfermions; in case of t, contribntions +m^ to the diagonal entries of the mass matrix 
also have to be added. Onr choice /x = 1 TeV together with the assnmption of gangino mass 
nnification ensnres that the LSP is an almost pnre bino. 

In contrast, in the second scenario we took fi = —100 GeV, rrig = 800 GeV, so that the LSP 
is dominated by its higgsino components, althongh the bino component still contribntes ~ 20%. 
In this scenario we took rriq = 1.5 TeV for all sqnarks and mj = 1.2 TeV for all sleptons, since 
in mSUGRA large scalar masses are reqnired if the LSP is to have a large higgsino component. 
We took GP-odd Higgs boson mass tua = 1 TeV in both cases, and tan/? = 3.6; we jnst saw 
that the latter choice is not important for ns. In the following we will refer to these two choices 
as the “gangino” and “higgsino” set of parameters, respectively. 

In Fig. ini we compare the FFs of an initial first or second generation SU{2) donblet qnark 
qi for these two scenarios. The main difference occurs in the FF into the LSP, which is 
significantly softer for the higgsino set. The reason is that most heavy superparticles (sfermions 
and gluinos) preferentially decay into gaugino-like charginos and neutralinos, which have much 
larger couplings to most sqnarks than the higgsino-like states do. These gaugino-like states 
are the lighter two neutralinos and lighter chargino in case of the gangino set, but they are the 
heavier x states for the higgsino set. The supersymmetric decay chains therefore tend to be 
longer for the higgsino set, which means that less energy goes into the LSP produced at the 
very end of each chain. 

Fig. U\ shows the same comparison for an initial first or second generation SU (2) doublet 
squark q^. Not surprisingly, the FFs of a squark are more sensitive to details of the sparticle 
spectrum than those of a quark. In particular, in addition to the reduced FF into the LSP, 
we now also see that the FFs into neutrinos and electrons are suppressed for the higgsino set 
relative to the gangino set. This is partly again due to the longer decay chains, which pushes 
these FFs towards smaller x where the normalization factor suppresses them more strongly, 
and partly because the branching ratios for leptonic decays of the SU (2) gaugino-like x states 
are smaller here than for the gangino set, which implies that fewer leptons are produced in 
sparticle decays. On the other hand, the longer decay chains and larger hadronic branching 
ratios for x decays characteristic of the higgsino set lead to an increase of the total multiplicity 
of 25% or so, as can be seen from the FFs at small x; of course, in this region the ratios of 
these FFs again approach their universal values, as discussed in Sec. 3.1. 

Hn our case X particles could contribute significantly to the Dark Matter; in this scenario, which is realized 
only for a small region of the total allowed Mx , Tx plane, the upper bound on the LSP relic density would have 
to be tightened accordingly, but the allowed regions of parameter space would be qualitatively the same. 
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If the initial particle is strongly interacting, the rapid evolution of the shower ensures that 
the generalized FFs m describing the evolution between Msusy and Mx essentially vanish 
at X ~ 1, i.e. all spectra are smooth. In contrast, if the initial particle I has only weak 
interactions, a signihcant 5—peak will remain at x = 1 in the generalized FF Dj. If / is a 
superparticle or Higgs boson, the decays of I can therefore lead to sharp edges in the hnal FFs. 
This is illustrated in Fig. |H1 which shows the FFs for an initial hrst or second generation SU{2) 
doublet sleptons II- The parameters of the gaugino set are chosen such that II sleptons can 
only decay into I + LSP. The decays of the II which survive at x = 1 therefore lead to edges 
in the FFs into e, Ue and recall that II is an equal mixture of e^, 9^, jli and z/^. The edge 
in the FF into e occurs at a somewhat larger value of x than those in the FFs into since 
after SU{2) symmetry breaking the charged members of the slepton doublets are a little heavier 
than the neutral ones; the decay cl —^ ex? therefore deposits more energy in the electron than 
z/g —>■ z/gXi deposits in the neutrino. However, in both cases the bulk of the energy goes into the 
LSP, which is rather close in mass to the slepton. This is quite different for the higgsino set, 
where sleptons are much heavier than all x states. As a result, almost the entire slepton energy 
can go into the decay lepton, leading to FFs into e, z/g and z/^ that are peaked very near x = 1 
(after multiplying with x^). Furthermore, since most sleptons now hrst decay into heavier x 
states rather than directly into x?; the FF into the LSP is much softer than for the gaugino set. 
Finally, the effect of the longer decay chains of SUSY particles on the overall multiplicity now 
amounts to about a factor of 2, and is thus much more pronounced than for initial squarks; 
this can be explained by the reduced importance of the shower evolution in case of only weakly 
interacting primaries. 

Fig. ini shows that in case of an initial Hi Higgs doublet, the role of the two parameter sets 
is in some sense reversed. Recall that we chose tan/? > 1 and 3> Mz- In that case the 
heavy Higgs bosons mostly consist of various components of the Hi doublet, with only small 
admixtures of H 2 ] see eq. JED in Appendix B. As usual with only weakly interacting primaries, 
the generalized FF remains sizable at x = 1 even at scale Mx- In the higgsino set, the 
dominant decay modes of the heavy Higgs bosons involve a gaugino and a higgsino, leading to 
a large FF into the LSP in this case. Since in the gaugino set the mass of the higgsino-like x 
states is very close to the mass of the heavy Higgs bosons, these supersymmetric decay modes 
are closed for the heavy Higgs bosons in this case, which instead predominantly decay into top 
quarks, with decays into b quarks and r leptons also playing some role. The fragmentation 
and decay products of these heavy quarks lead to a signihcantly larger FF into protons in the 
gaugino region; semi-leptonic t and b decays as well as the r decays also lead to enhanced FFs 
into electrons and neutrinos for the gaugino set. Finally, the hadronic showers initiated by the 
decay products of the top quarks as well as by the b quarks produced directly in the decays of 
Higgs bosons raise the total multiplicity for the gaugino set to a value which is slightly larger 
than that for the higgsino set. 

As hnal example we compare the FFs of an initial H 2 higgsino doublet in Fig. Here we 
again hnd a larger FF into the LSP for the higgsino set, including a peak at x = 1. In this case 
this is simply a rehection of the large iL? component of the LSP. On the other hand, in case of 
the gaugino set H 2 projects almost exclusively into the heavier x states, which have many two- 
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body decay modes into sleptons and leptons. This explains the relative enhancement at large 
X of the FFs into leptons that we observe for the gaugino set, as well as the structures in these 
FFs. On the other hand, the longer sparticle decay chains again imply a somewhat larger 
overall multiplicity for the higgsino set. These decays of heavy sparticles are important here 
since the large top Yukawa coupling of H 2 initiates a signihcant parton shower in this case, 
where numerous superparticles are produced. This is quite different for an initial Hi at small 
tan/3 (not shown), where we hnd a sma//er overall multiplicity for the higgsino set, since the 
number of produced superparticles remains small, and the initial particle Hi has a longer decay 
chain for the gaugino set. 

Altogether we see that the SUSY spectrum can change the hnal FFs, and thus the hnal 
spectra of X decay products, signihcantly. Generally this effect is stronger for an initial super¬ 
particle or heavy Higgs boson than for an SM particle, and stronger for only weakly interacting 
particles than for those with strong interactions. However, with the exception of the FFs into 
the LSP, the variation is usually not more than a factor of two, and often much less. The 
dependence of the X decay spectra on SUSY parameters can therefore be signihcant for de¬ 
tailed quantitative analyses, but this dependence is always weaker than the dependence on the 
primary X decay mode(s). 

3.4 Coherence effects at small x: the MLLA solution 

So far we have used a simple power law extrapolation of the hadronic (non-perturbative) FFs 
at small x. This was necessary since the original input FFs of ref. Ea are valid only for x > 0.1. 
As noted earlier, we expect our treatment to give a reasonable description at least for a range 
of X below 0.1. However, at very small x, color coherence effects should become important |02j . 
These lead to a battening of the FFs, giving a plateau in xD{x) at Xpiateau ~ \jQhiid/Mx ~ 10“® 
for Mx = 10^® GeV. One occasionally needs the FFs at such very small x. For example, 
the neutrino hux from X decays begins to dominate the atmospheric neutrino background at 
E ~ 10^ GeV I2ZIEH1, corresponding to x ~ 10“^^ for our standard choice Mx ~ 10^® GeV. In 
this subsection we therefore describe a simple method to model color coherence effects in our 
FFs. 

This is done with the help of the so-called limiting spectrum derived in the modihed leading 
log approximation. The key difference to the usual leading log approximation described by the 
DGLAP equations is that QGD branching processes are ordered not towards smaller virtualities 
of the particles in the shower, but towards smaller emission angles of the emitted gluons; note 
that gluon radiation off gluons is the by far most common radiation process in a QGD shower. 
This angular ordering is due to color coherence, which in the conventional scheme begins to 
make itself felt only in NLO (where the emission of two gluons in one step is treated explicitly). 
It changes the kinematics of the parton shower signihcantly. In particular, the requirement that 
emitted gluons still have sufficient energy to form hadrons strongly ahects the FFs at small x. 
For sufficiently high initial shower scale and sufficiently small x the MLLA evolution equations 
can be solved explicitly in terms of a one-dimensional integral [22] • This essentially yields the 
modihed FF describing the perturbative gluon to gluon fragmentation, in the language of 
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eq. (uni). In order to make contact with experiment, one makes the additional assumption that 
the FFs into hadrons coincide with D^, up to an unknown constant; this goes under the name 
of “local parton-hadron duality” (LPHD) jSSI- Here we use the £t of this “limiting spectrum” 
in terms of a distorted Gaussian [21] , which (curiously enough) seems to describe LEP data on 
hadronic FFs somewhat better than the “exact” MLLA prediction does. It is given by 
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where rij is the average multiplicity. The other quantities appearing in eq. m are dehned as 
follows: 
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where (3 is the coefEcient in the one-loop beta-function of QCD and p = 11 + 2Nf/27, Nf 
being the number of active flavors. Eqs. m and (P|) have been derived in the SM, where 
/3 = 11 - 2Nf/3. Following rehUH we assume that it remains valid in the MSSM, with /3 = 3 
above the SUSY threshold Msusy and p = 11 + 8/9. Note that we do not attempt to model 
the transition from the full MSSM to standard QCD here; indeed, we do not know of an easy 
way to do this, since the limiting spectrum cannot be written as a convolution of two other 
spectra. On the other hand, the position ^ of the plateau depends only on and only via 
the second term, which is suppressed by a factor yT ~ 6.5, whereas the parameters a and s 
describing the behavior in the vicinity of the maximum depend in leading order in r only on 
Finally, the coefEcient p is very similar in the SM and MSSM. We therefore expect the 
error we make by ignoring the transition from MSSM to SM to be smaller than the inherent 
accuracy of eq. eg. 

When comparing MLLA predictions with experiments, the overall normalization Uj (which 
depends on energy) is usually taken from data. We cannot follow this approach here, since 
no data with Q ~ Mx are available. Moreover, usually MLLA predictions are compared with 
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Figure 11: Comparison between the MLLA solution and our results without coherence effects, 
for the final proton and neutrino spectra. We assume that X undergoes two-body decay into 
(IlQi- 

inclusive spectra of all (charged) particles. We need separate predictions for various kinds 
of hadrons, and are therefore forced to make the assumption that all these FFs have the 
same x—dependence at small x. This is perhaps not so unreasonable; we saw above that the 
DGLAP evolution predicts such a universal x—dependence at small x. We then match these 
analytic solutions (HHD, (HI with the hadronic FFs we obtained from DGLAP evolution 
and our input FFs at values Xq, where for each hadron species h the matching point Xq and the 
normalization fih are chosen such that the FF and its first derivative are continuous; we typically 
find Xo ~ 10“^. Note that this matching no longer allows to respect energy conservation exactly. 
However, since the MLLA solution begins to deviate from the original FFs only at a; ~ 10“^, 
the additional “energy losses” are negligible. 

Some results of our MLLA treatment are shown in Fig. [TTJ Here the “non-MLLA” curves 
have been obtained by extrapolating our numerical results described earlier, which extend 
“only” to a; = 10“^, by using simple power-law fits. We see that at a; ~ 10“^^ the FFs are 
suppressed by about two orders of magnitude, but the effect diminishes quickly at larger values 
of X. Note that the FFs into protons and into neutrinos have slightly different shapes in the 
small—X region. By assumption the FFs have the same shape for all hadrons; however, in going 
from the spectrum of pions and kaons to the neutrino spectrum, several additional convolutions 
are required, which shift the peak of the distribution to even smaller values of x. This figure also 
shows that the MLLA predictions closely tracks the non-MLLA solution for x values that are 
several orders of magnitude smaller than the matching point xq; this illustrates the advantage 
of requiring both the FF and its first derivative to be continuous at xq. 
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4 Summary and Conclusions 

In this article, we presented a detailed analysis of the decay of a very massive X particle, 
extending onr earlier work ng. In particular, we were able to improve the accuracy of our 
code; as a result, we can now ensure energy conservation to a numerical accuracy of better 
than 1%, as compared to up to several % in ref. ng. Moreover, we showed that the dependence 
of our results on the necessary extrapolation of the measured fragmentation functions (FFs) 
towards small x is negligible. We also included leading higher-order QCD corrections at very 
small X using the MLLA approximation for taking into account color coherence effects; this 
approximation is in good agreement with data from particle colliders. These effects become 
significant for x ^ 10“^, decreasing the predicted fluxes at a: ~ 10“^^ by about two orders of 
magnitude. 

Furthermore, we showed that varying SUSY parameters can have some impact on our 
results, affecting the shapes of the FFs at a: > 0.01 and in some cases also the total multiplicity; 
however, the dependence on the SUSY spectrum is much milder than the dependence on the 
primary X decay mode(s). Qualitatively the photon and LSP fluxes are the most important 
ones at large x if the primary is a strongly interacting (s)particle; if the primary has only weak 
interactions, the lepton fluxes can also be very large at large x. The proton flux is always 
subdominant in this region. In contrast, the shapes of most FFs at small x can be predicted 
almost uniquely. This leads to the following ordering of the fluxes at x < 0.01: the largest flux 
is of muon neutrinos, followed by photons, Ug and electrons, and finally protons. The ratios of 
these fluxes become almost independent of x in this region, the proton flux being about a factor 
of five smaller than the flux. On the other hand, the two smallest fluxes at small x, of LSPs 
and finally do depend sensitively on various currently unknown parameters. Generically 
they rise less rapidly with decreasing x than the other fluxes do; already at x 10 the Ut 
and LSP flux are usually about one order of magnitude below the proton flux. 

Finally, in the appendices we give additional details of our description of the complete 
cascade. In particular. Appendix A contains the first complete set of leading order splitting 
functions for the MSSM, including all gauge as well as third generation Yukawa interactions. 
A “catalogue” containing an almost complete set of FFs for a given set of parameters is given 
in Appendix E. 

This work presents the to date most accurate and complete description of the spectra at 
source of stable particles resulting from the decay of a superheavy X particle. These spectra 
are needed for all quantitative tests of the “top-down” explanation of the most energetic cosmic 
ray events. Of course, in order to be able to compare with fluxes measured on or near Earth, 
effects due to the propagation through the galactic, and perhaps extragalactic, medium jH] have 
to be included, which depend on the distribution of X particles throughout the Universe; we 
have made no attempt to do this. On the other hand, our description of X decays is model- 
independent in the sense that it allows to incorporate any primary X decay mode. Indeed, it 
could with very little modification also be used to describe the evolution of very energetic jets 
produced through some other mechanism (e.g. the annihilation of very massive stable particles), 
as long as the initial virtuality of the produced particles is comparable to their energy. 
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Returning to the original problem of ultra-high energy cosmic rays (UHECR), the biggest 
obstacle towards a test of generic top-down models is the strong dependence of the predicted 
decay spectra on the primary decay mode. Most previous investigations assumed that X decays 
into a pair of quarks, but we are not aware of any compelling argument why this should be 
the dominant decay mode. On the other hand, data may already rule out some classes of 
top-down models. For example, it seems likely that few, if any, UHECR are photons jH]. In 
the context of top-down models, this leaves protons as only choice. Our results then seem 
to disfavor models where X decays primarily into particles with only weak interactions, since 
this implies a large ratio of the photon to proton flux at large x. However, this argument 
may not apply if Mx ^ 10^^ GeV, since then all events seen so far are aX x % 0.01, where 
the ratio of photon to proton fluxes is essentially independent of the primary X decay modes. 
Moreover, the photon flux may be diminished more efficiently between source and detector 
than the proton flux. Searches for very energetic neutrinos might therefore lead to somewhat 
more robust tests of top-down models 1211 EH; as noted earlier, the predicted neutrino flux 
should begin to exceed the background from atmospheric neutrinos at very small values of x. 
Nevertheless, the need to normalize the expected flux to the observed flux of UHECR events, 
and hence to the proton and perhaps photon flux at much larger x, re-introduces a large model 
dependence even in this case |2S]- Moreover, other proposed explanations of the UHECR also 
predict sizable neutrino fluxes at very high energy, e.g. due to the GZK process itself. The 
failure to observe such neutrinos could therefore exclude top-down models (given sufficiently 
large detectors), but a positive signal may not be sufficient to distinguish them from generic 
“bottom-up” models. This discrimination might be achieved by searching for the predicted 
flux of very energetic LSPs, since the LSP flux in bottom-up models is undetectably small; 
however, this test will require very large detectors |121- We conclude that ultimately the test 
of this idea will probably require a combined analysis of different signals, at quite different 
energies and in different detectors. We provide one of the tools needed to perform such an 
analysis, since we are able to systematically study the fluxes of all stable particles at source, 
and their correlations, for all top-down models. 
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Appendix A: Splitting functions of the MSSM 

The splitting function (SF) Pji{x) describes the radiation of particle j off particle i. Its 
X—dependence is determined by the Lorentz structure of the corresponding vertex, while the 
normalization also depends on the associated group [color and SU{2)] factors. If there is no 
vertex relating these two particles the SF is simply 0. We hrst list the functional forms we will 
need, together with the spins of the particles involved in the branching process i ^ j + k {V 
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for vector, F for spin—1/2 fermion, S for scalar): 
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For convenience, we also define (1') = (1) + (0) and (7') = (7) + (0). 
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Table 4: SUSY-QCD splitting fnnctions Pji, where j and i determine the row and colnmn 
of the table, respectively. The fnnctional forms denoted by (n), n = 0,..., 10 have been 
dehned in eq. (HU), with (1') = (1) + (0) and (7') = (7) + (0). The “mnltiplicity factors” are: 

= 1, Ntj^ = Nuj^ = Ndji = 2 and = 4. In the MSSM phase, i.e. for Q > Msusy, 
the nnmber of active flavors (quarks and squarks) is F = 6. 


The 16 SFs of SUSY-QCD listed in table IH are derived from [22]; in eq. (fn|l they come 
with a factor of the strong coupling as- Note that in ref. I2ZI the chirality index L, R was 
always summed over; e.g. Ptg = Pt^^g + Pt^g, where ti now only describes the left-handed 
top quark (and not the third generation quark doublet). Since these two terms are equal. 
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one has Pt^^g = Ptj^g = Pqg/‘^- On the other hand, our “(s)quark” distributions always include 
anti(s)quarks. This re-introduces a factor of 2, so that for us e.g. Ptj^g = Pqg of j2Zj. Additional 
factors arise for (s)quarks of the hrst and second generation. As described in Sec. 2.2, we always 
average over (s)quarks and anti(s)quarks with given hypercharge of the hrst two generations. 
This implies Pu^g = Pd^g = ‘^Pqg and Pq^g = ^Pqg-, where the additional factor of two in the 
second expression comes from summing over the SU{2) index of the doublet qi. The same 
factors appear in SFs describing gluon to squark splitting as well as gluino splitting into a 
squark and a quark. A complete list of these factors Nq is given in the table caption. On the 
other hand, in the absence of havor-changing interactions SFs involving quarks and squarks 
only always come with factor 1 if the “compound particles” ur, ql etc. are properly normalized. 

The SFs stemming from electroweak interactions have similar structures; we just need to 
compute the correct group and multiplicity factors. The results are listed in tables El and IHl In 
these tables we list SFs including the appropriate multiplicity factor; a single SU{2) doublet 
without antiparticles would have Nf = 1/2. Note that there is no difference between Higgs and 
SU{2) doublet lepton superhelds as far as gauge interactions are concerned. Finally, due to 
the absence of self-interactions of 17(1) gauge bosons, the splitting functions Pbb and are 
pure delta-functions, with coefficients fixed by energy conservation, eq. ©. In all three gauge 
interactions we find that the coefficient of the <5—function is the same in Pff and ^ff for any 
matter fermion /, and also in Pyy and Pvv for a gauge boson V; this latter coefficient is —1/2 
times the coefficient in the /3—function of the corresponding gauge coupling. 
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Table 5: SU{2) splitting functions Pji, where particles j and i are associated with the row and 
column, respectively. The functional forms denoted by (n), n = 0,..., 10 have been defined in 
ea. IjA.Hl . with (1') = (1) + (0) and (7') = (7) + (0). Nd is the total number of SU{2) doublets; 
in the MSSM, Nd = 14. / stands for any matter or Higgs fermion, with Nf being the number 
of doublets (not counting anti-doublets) described by fi or //,. For our “compound” states, 
these are: Nq^ = 6, Ni^ = 2, = 3, N^^ = Nh^ = Nh^ = 1. 


Finally, Yukawa couplings only appear in P/l/r, hfifji and hfifn vertices. We therefore 
only need functional forms (3), (4) and (5) from ea. IlA.H) . The coefficients can be determined 
from the analogous terms due to U{1)y interactions by replacing (gyYfY by Aj/2, where the 
extra factor of 1/2 corrects for the factor ^/2 appearing in front of gaugino-fermion-sfermion 
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Table 6: U{1)y splitting functions Pjt, where particles j and i are associated with the row and 
column, respectively. The functional forms denoted by {n), n = 0,..., 10 have been dehned in 
eq. (D, with (1') = (1) + (0) and (7') = (7) + (0). The sum of squared hypercharges of all 
particles S/ = H iii fhe MSSM. / stands for any matter or Higgs fermion with hypercharge 
Yf, while Uf is the number of degrees of freedom (not counting anti-particles) described by / 
or /. For our “compound” states, these are: = 1/36, Uq^ = 12; Y//^ = 4/9, = 6; Y^^ = 

l/9,nrf^ = 6; = 1/4, = 4; Y/i = l,ni^ = 2; = 1/36, = 6; = 4/9, = 3; 


Ybn = ^bn = 3; F4 = F^^ = F^ = 1/4, = uh, = rin^ = 2; F^ = 1, = 1. 
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Table 7: Splitting functions Pji originating from Yukawa interactions, where particles j and 
i are associated with the row and column, respectively. The functional forms denoted by 
(n), n = 0,3,4, 5 have been dehned in eq. (HH). Since we only include Yukawa interactions 
for the third generation, we only have to consider three cases. For the top Yukawa coupling, 
fi = tL, Jr = tji, H = H2 and number of colors Nc = 3; for the bottom Yukawa coupling, 
= ti) Ir = bR, H = Hi and Nc = 3; hnally, for the tau Yukawa coupling, fL = TL-ifR = 
tr, H = Hi and Nc = 1. 


37 


























vertices in the supersymmetric Lagrangian. Since Yukawa interactions couple matter fields 
with different chiral indices all diagonal SFs due to Yukawa couplings are pure 5—functions, the 
coefficients again being determined by energy conservation; as before, we find equal coefficients 
for diagonal SFs of a particle and its superpartner. The resulting SFs are listed in table [Tj 
As usual, these SFs are multiplied with Q!//(27r) = Aj/(87r^) in the DGLAP equations. The 
three interactions we consider, involving the top, bottom and tan Yukawa couplings, can all be 
treated using table [3 by identifying the matter and Higgs fields appropriately and using the 
correct color factors, as explained in the caption. The SU (2) factors, which lead to the factor 
of 2 difference between SFs describing radiation off SU{2) singlet or doublet (s)fermions, are 
the same in all three cases. 1 

Appendix B: Unitary transformations between current 
and mass eigenstates in the MSSM 

In this Appendix we describe the unitary transformations occurring during the SUSY and 
SU (2) (g) U (1) breaking, where the quarks, leptons, weak gauge bosons and Higgs bosons as well 
as all superparticles acquire their masses jlHj • The superscript b denotes the mass eigenstates of 
the broken theory. The fields in the unbroken theory are the same as those described in Sec. 2.2. 
For example, stands for all left-handed quarks and antiquarks of the two first generations, 
i.e. the SU{2) doublets {uiyd^), {cl,sl) and their antiparticles (uiydL), (cl,sl), and thus 
describes eight degrees of freedom (times three, if color is counted separately). Similarly, 
II stands for both SU{2) doublets and and their antiparticles (e^,, i^e) and 

On the other hand, only describes u—quarks and their antiparticles, but includes 
both chirality states, and thus describes four degrees of freedom (not counting color). Recall 
that the transformation between mass and current eigenstates in eq. (USD only affects the upper 
index of the (generalized) FFs. In the given context qi therefore stands for the sum, not the 
average, of its “constituent fields”, as discussed in Sec. 2.2. Recall finally that massive gauge 
bosons “eat” Goldstone modes from the Higgs sector. These considerations lead to the following 
transformations for SM fields and Higgs bosons: 

b fell 

U = C = -qL + -Ur, 

d” = s‘ = i9L + io!H. 

b’^ = + bR, 

= -tr + tR, 


^Strictly speaking, Hi can only split into tr and tr, not into tr and tr, while the antiparticle H* can only 
split into tr and tr] analogous remarks hold for the other Yukawa-induced branching processes. However, this 
distinction plays no role for us, since we always average or sum over particle and antiparticle. 
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^L/i? = gL//? forg ='u,d,s,c, 

= \ (^b) h + sin^ (6»b) bn , 

^1 = ^cos^{9t)iL + sm‘^{0t)iR, 

62 = ^ sin^ {9b) ii + cos^ {9b) Ir , 

^ ^ sin^ {9t) ii + cos^ {9t) Ir , 

el = A = \ h , 

~ ^ 1 

— f^R — 2^^^ 

fi = ^ COs‘^{9r)fL + Sm^{9r)fR, 

fl = ^ sm^{9r) tl + cos^(6'^) tr , 
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= 




= 


9 = 


2 




Xl 


■= Xl + Xl = 


X2 •” X2 + X2 ~ 


Xi^ = 


sm^(jn) + sin2(7L)] y + cos^(7r) y + cos2(7l) y , 
cos2(7r) + cos^(7l)] y + sin2(7K) ^ + sin2(7L) y , 


B. 


( 2 ) 

2 Hi 


2H2 


2 W 

(i) 



IP 

^2 


7 ? ' 


V4 


(B.2) 


Here we have largely followed the notation of ISASUSY pg. However, we havejrsed the nwre 
common symbol x for charginos and nentralinos; in ISASUSY notation, Xi = fU_, X 2 — ^+1 
and xT = Zi- The mixing angles a (in the Higgs sector), 9i„ 6^, (in the sfermion sector), 
1li 1r (in the chargino sector) as well as the (in the nentralino sector) have been compnted 
numerically using ISASUSY. 


Appendix C: Two— and three—body decay spectra 


Generalities 

We want to define the decay functions (DFs) Pgs describing two- or three-body decay S —>■ s, 
see eq. (USD- These DFs can be obtained directly from the decay spectrum of S in the ultra- 
relativistic limit, where the energy Eg is much larger than the mass M of S': 


Pssiz) 


1 dTjEs) 
F dz 


(C.l) 


where z = Eg/Es- This spectrum can e.g. be evaluated by first computing the double dif¬ 
ferential decay distribution d?T/{dE^d cos 9*) in the rest frame of S', then boosting the four- 
momentum of s with boost factor 7 = Es/M at angle 9* relative to Ps, and finally integrating 
over cos 6** subject to the constraint that the boosted energy of s equals Eg. Note that eq. CD 
implies Jq Pgs(z)dz = 1; if S'—decays produce N identical particles s, the corresponding Pgs 
would thus have to be multiplied with an extra factor of N, in order to correctly reproduce the 
total multiplicity in the hnal state. Finally, momentum conservation implies Jfg Jq zPss{z) = 1. 

In case of two-body decays S —>■ i+j the energy E* in the rest frame of S' is hxed completely 
by the kinematics. The boost and integration over cos 9* then leads to a flat decay function: 



/mi+m2\‘^ 

V M J 



mi — m 2 

M 


2 


Q{z - z^f>) 0(4^ - z) (C.2) 
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for the decay product i with i = 1 or 2. The kinematic minimum and maximum of z are 
given by: 



mi + m2\ ^ 

M ) 


mi — m2 

M 


(C.3) 


For example, for mi —> M, m 2 —> 0, eq. (O implies —»• 1, z±^ —> 0, i.e. the entire energy 

of S goes into the massive decay product. In contrast, for mi = m 2 , the energy of S will on 
average be shared equally between the two decay products; if mi = m 2 —> 0, the 2 ;^*^ can lie 
anywhere between zero and one. Since E* is hxed in this case, our treatment of two-body 
decays is exact up to possible polarization effects; we do not expect these effects to be very 
important, except perhaps in case of r decays (which, however, usually do not contribute very 
much to the hnal spectra of stable particles). 

Three-body decays lead to a nontrivial distribution of the energy of the decay products 
already in the rest frame of S. For simplicity we assume that at most one of the three decay 
products is massive; this should be a safe approximation, except for b ctUt decays, which 
have a rather small branching ratio. We then need separate DFs for the massive and massless 
decay products. For the massive decay product, with mass m, we hnd 




. 1 


(C.4) 


where 2 ; G [^, 1] and the normalization factor is given by: 


N. = 


1 / \ 

2V 


-1 


(C.5) 


If on the contrary, s is one of the massless decay products, we find: 



N 3 



9 

m 

jp 


1 + log — + log(l 


m^ 



(C.6) 


where now z G [0,1 — ^]; the normalization factor N 3 has been given in ea. ()(h5fl . 

Our treatment of three-body decays is not exact, since it ignores dynamical effects (de¬ 
scribed by the invariant Feynman amplitude) on the decay spectrum in the S rest frame.* 
However, treating these effects properly is quite nontrivial, since it would force us to introduce 
many different three-body decay functions. Note in particular that massive superparticles 
(charginos and neutralinos) do generally not decay via V — A interactions, unlike the b and c 
quarks and heavy /i and r leptons in the SM. Moreover, the Feynman amplitudes in many cases 
depend nontrivially on the polarization of the decaying particle; this could only be described 

*The calculation of the corresponding branching ratio in ISASUSY does include these dynamical effects; in 
other cases the required branching ratio can be taken directly from experiment, e.g. for r decays. 
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at the cost of introducing many additional generalized fragmentation functions, since we would 
have to keep track of left- and right-handed particles separately. However, experience from 
hadron collider physics teaches us that including the exact decay matrix elements is usually not 
very important if one is only interested in single-particle inclusive spectra. We expect this to 
be true in our case as well, since the convolution with parton distribution functions necessary at 
hadron colliders is reminiscent of the convolution with generalized FFs in our case. We hnally 
note that longer decay chains involving two- and three-body decays can be treated by simply 
convoluting appropriate factors of and ■ 

Treatment of heavy quark decays 

The top quark being very heavy {mt ~ 175 GeV 3> ruhad ~ 1 GeV), it decays before hadroniz- 
ing, and can thus be included in the decay cascade at scale Msusy- On the other hand, the 
hadronization of the b and c quarks has to be treated with some care. The “input” fragmen¬ 
tation functions we used (22] already include the hnal hadrons (nucleons, kaons and pions) 
produced at the end of the decay cascade of c— and 6—flavored hadrons. However, they do 
not include the leptons arising from this cascade, which are not negligible. We therefore imple¬ 
mented a special treatment for this component, using the empirical FFs proposed by Peterson 
et ah [22] for heavy quarks as a basis for the fragmentation of c— and 6—hadrons. To that 
end, we used two “generic” particles, a c— and a 6—hadron, with respective average masses 
rUc = 2.1 GeV and mb = 5.3 GeV; we also had to renormalize the complete set of FFs for 6’s 
and c’s. The scheme can be described by Fig. 12. Here, Bi{b) and Bi{c) describe the branching 
ratio of the semi-leptonic decay modes of b— and c—flavored hadrons, respectively [summed 
over all accessible pairs 


b -*- 5—hadron 

Peterson 


semi-leptonic 
decay P/(5) 


c—hadron I i>i 

semi-leptonic 

-►--hadron I ui 

decay P/(c) 


c -- c—hadron 

Peterson 


semi-leptonic 
decay Bi{c) 


s—hadron I Ui 


Figure 12: Schematic hadronization and decay cascade for heavy quarks c and b. The 
“s—hadrons”, mainly kaons, are already included in the FFs given in |25j . 
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As mentioned earlier, the leptonic b and c decay products have to be included in the nor¬ 
malization of the FFs and To that end, we introduce Rc and Rg, the energy carried by 
the c and s quark in semi-leptonic b— and c— decays, respectively, as well a.s xb and xd, the 
energy fraction of the b (c) quark carried by the 6—flavored (c—flavored) hadron. The latter 
are given by 


xb,d 





Pet 


[z) dz , 


(C.7) 


where is the Peterson FF [321 for b and c quarks, respectively; we took Cc = 0.15, et = 
0.015 for the single free parameter in these FFs. We compute Rc and Rg from pure phase 
space, i.e. we again ignore possible effects of the Feynman amplitudes on the three-body decay 
distributions. This gives: 


X ^ dTib^du) 

mfeF(6 —t du) Jmc ^ dEc 

(1 - rf 

3(1 — r^) — 6r logr ’ 


(C.8) 


where r = = 0.157 for our choice of average b— and c—hadron masses; note that 

Rc —t 1/3 (1) for r —t 0 (1). Ea. (l(h8fl can also be used for the computation of Rg, with 
uTc —t rUg ~ 0.5 GeV, rrib —t We- The FFs of jlH] only include the hadrons produced in the 
fragmentation and decays of the c— and 6—quarks. Their normalization, which we need to 
know for the necessary extrapolation of these FFs towards small x as discussed in Appendix 
D, is therefore given by 


1 

Xb 


E 



zD^{z) dz 


+ 

+ 

+ 

+ 


I- Xb 
Xb 

[l-Bi{b)]-[l-Bi{c)] 

Bi{b)-[l-Bi{c)]-Rc 

[1 ~ Bi{h)\ ■ Bi(c) ■ {Rc ■ Rg 1 — Rc) 

Bi{b) ■ Bi{c) ■ Rc- Rg 

E _ Bi{b) • (1 - Rc) - Rc ■ Bi{c) ■ (1 - Rg). 


(C.9) 


The right hand side can be understood as follows: the hrst line describes the contribution of the 
light hadrons produced when the 6—quark hadronizes into a 6—flavored hadron; the second line 
describes purely hadronic decays; the third line describes leptonic primary b decays followed by 
hadronic c decays (in this case only the fraction Rc of the 5—hadrons energy goes into hadrons); 
the fourth line describes the hadronic energy fraction in the case of a hadronic primary b decay 
followed by leptonic c—decays; hnally, the fifth line describes the hadronic energy fraction after 
leptonic decays in both the primary and secondary decays. The same holds for c—hadron 
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decays, up to the simplifying fact that “s—hadrons” are already fully included into the FFs of 
|2ni- We get: 

— ^ zD^^iz) dz = — - Bi{c) + Bi{c) ■ R,. (C.IO) 

Xd ^ Jo Xd 

Appendix D: Parameterization of the input fragmentation 
functions 

Here we give the input fragmentation functions (FFs) we used to describe the hadronization of 
quarks and gluons, taken from |2SI, and the parameters of the extrapolation we made in the 
small X region. 

The functions taken from ref [23] are given with the functional form — x)^ ; there are 

given in table |Hl at the scale where quarks and hadrons hadronize, i.e. Qq = max{mq,Qhad)- 
We used the NLO results, excepted for the s quark, for which we had to use the LO ones, 
because the NLO form didn’t allow us to impose energy conservation and continuity at low x. 

As we showed in Sec. ESI the hnal result at low x depends very little on the chosen power 
law in our parameterization 

f{x) = ax~°‘ + b\ogx + c, a > 0, (D.l) 

once energy conservation has been imposed. Here we therefore only give results for a pa¬ 
rameterization where a' is taken to be 1. That is, we assume that the multiplicity due to 
non-perturbative effects gets the same contribution for each decade of energy, if the hadron’s 
energy is small compared to that of the initial parton. 

In order to obtain a unique solution with the only two constraints at our disposal (energy 
conservation and continuity of the FFs), we imposed these constraints on the sum of the FFs 
where z is a given initial parton, and h runs over the hnal hadrons. Of course, 
energy will be conserved independently for each initial parton i. For each i we dehne a cut¬ 
off Xq which dehnes the transition between the functions given in |23| and our extrapolation. 
The Xq have to be chosen such that the equations of energy conservation admit a solution; 
a necessary (but generally not sufficient) requirement is that the integral over the original 
FFs satisfy Jl^dzz D^{z) < 1. Our requirement of continuity at Xq implies that the hnal 
results depend very little on the precise values of the Xq. For simplicity we assume that all the 
D^{x, Qhad) ^ given i have the same shape at small x] recall that purely perturbative ehects 
ensure that this is true after DGLAP evolution, which anyway greatly reduces the sensitivity 
to the input. The normalizations for the various hadrons can then be read oh directly from the 
results of ref. once Xq has been determined. The results are presented in table. El which 
lists the cut-oh Xq, the coefficients a, b, c of the functional form EH describing the sum J2h 
for hxed parton i, and the normalization coefficients^ iV/*, so that of course, 

EhNt = iyt. 


tAt first sight the relative ordering of the , Nf 


factors may seem counter-intuitive. Indeed, a u—quark 
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Dpix,Ql) 

= Nx°'{l — xY 

P 

n 

TT^ 

7r° 



u 

N = 1.26 

0.63 

0.448 

0.224 

0.178 

4.96 


a = 0.0712 

0.0712 

-1.48 

-1.48 

-0.537 

0.0556 


[3 = 4.13 

4.13 

0.913 

0.913 

0.759 

2.8 

d 

0.63 

1.26 

0.448 

0.224 

4.96 

0.178 


0.0712 

0.0712 

-1.48 

-1.48 

0.0556 

-0.537 


4.13 

4.13 

0.913 

0.913 

2.8 

0.759 

s 

4.08 

4.08 

22.3 

11.15 

0.259 

0.259 


-0.0974 

-0.0974 

0.127 

0.127 

-0.619 

-0.619 


4.99 

4.99 

6.14 

6.14 

0.859 

0.859 

c 

0.0825 

0.0825 

6.17 

3.085 

4.26 

4.26 


-1.61 

-1.61 

-0.536 

-0.536 

-0.241 

-0.241 


2.01 

2.01 

5.6 

5.6 

4.21 

4.21 

b 

24.3 

24.3 

0.259 

0.1295 

1.32 

1.32 


0.579 

0.579 

-1.99 

-1.99 

-0.884 

-0.884 


12.1 

12.1 

3.53 

3.53 

6.15 

6.15 

9 

1.56 

1.56 

3.73 

1.865 

0.231 

0.231 


0.0157 

0.0157 

-0.742 

-0.742 

-1.36 

-1.36 


3.58 

3.58 

2.33 

2.33 

1.8 

1.8 


Table 8: Input fragmentation functions at small x, with functional form iVx“(l — x)^, taken 
from [23] at Qq = max{mq,Qhad)- We took their NLO results for u,d,c,b and g, but the LO 
result for the s quark. See the text for further details. 


initial parton 

Xo 

a 

b 

c 

P 

n 

TT^ 

7r° 



u 

0.27 

4.06 

-9.74 

-14.40 

0.05 

0.025 

0.38 

0.19 

0.05 

0.31 

d 

0.27 

4.06 

-9.74 

-14.40 

0.025 

0.05 

0.38 

0.19 

0.31 

0.05 

s 

0.20 

5.74 

-18.47 

-31.42 

0.14 

0.14 

0.41 

0.21 

0.05 

0.05 

c 

0.27 

4.06 

-4.16 

-6.24 

0.05 

0.05 

0.30 

0.15 

0.22 

0.22 

b 

0.20 

5.74 

-27.81 

-49.27 

0.08 

0.08 

0.35 

0.17 

0.16 

0.16 

9 

0.37 

1.82 

-4.81 

-2.40 

0.05 

0.05 

0.50 

0.25 

0.07 

0.07 


Table 9: Coefficients of the extrapolation of the input fragmentation functions at small x. 
Column 2 gives the cut-off Xq where we switch from the FFs from [23] to a parameterization 
in the form (jP.ljl . Columns 3 to 5 give the coefficients of this parameterization, as applied 
to the sum The remaining columns give normalizations so that 

Note that h always stands for the sum of particle and anti-particle, whenever the two are not 
identical; for example, vr^ stands for the sum of and 7r“, stands for the sum of and 
K^, etc. 
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Appendix E: Stable particle spectra for different initial 
(super)particles 

Here we give an almost complete set of FFs for different initial particles, for one set of SUSY 
parameters, with low tan/3 and gangino-like LSP; the dependence of these results on the SUSY 
parameters has been analyzed in Sec. Id.dl We used a ratio of Higgs vevs tan/3 = 3.6, a gluino 
and scalar mass scale Mgusv ~ 500 GeV, a supersymmetric Higgs mass parameter /i = 500 
GeV, a GP-odd Higgs boson mass uia = 500 GeV and trilinear soft breaking parameter A = 1 
TeV. As usual, we plot ■ Df {x, Mx)- We take Mx = 10^® GeV, as appropriate for a GUT 
interpretation of the X particle. 


should more readily fragment into a charged Kaon than into a neutral one, whereas the opposite behavior is 
expected for d—quarks. Recall, however, that here we are only interested in the behavior at small x. In this 
case ref. m finds the opposite behavior as at large x, i.e. rt—quarks indeed seem to be more likely to produce 
a soft neutral kaon than a soft charged kaon. 
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